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Abstract. In this article, we discuss the stability of soft quasicrystalline phases in a 
coupled-mode Swift-Hohenberg model for three-component systems, where the char¬ 
acteristic length scales are governed by the positive-definite gradient terms. Classic 
two-mode approximation method and direct numerical minimization are applied to 
the model. In the latter approach, we apply the projection method to deal with the po¬ 
tentially quasiperiodic ground states. A variable cell method of optimizing the shape 
and size of higher-dimensional periodic cell is developed to minimize the free energy 
with respect to the order parameters. Based on the developed numerical methods, 
we rediscover decagonal and dodecagonal quasicrystalline phases, and find diverse 
periodic phases and complex modulated phases. Furthermore, phase diagrams are 
obtained in various phase spaces by comparing the free energies of different candidate 
structures. It does show not only the important roles of system parameters, but also 
the effect of optimizing computational domain. In particular, the optimization of com¬ 
putational cell allows us to capture the ground states and phase behavior with higher 
fidelity. We also make some discussions on our results and show the potential of ap¬ 
plying our numerical methods to a larger class of mean-field free energy functionals. 
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1 Introduction 

Since the first discovery of quasicrystal in a rapidly-quenched Al-Mn alloy in 1980s lU, 
hundreds of metallic quasicrystals, whose building blocks are on the atomic scale, are 
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found both in syntheses in the laboratory ||2]-[3 and in nature ||8|I. More recently, a growing 
number of mesoscopic quasicrystalline orders are reported in soft-matter systems 19141211. 
These soft quasicrystals exhibit distinct properties from their solid-state counterparts. 
For example, dodecagonal or even higher order symmetry are possible in the soft-matter 
systems ll9l [T0l[T3ll , while dodecagonal solid-state quasicrystals are rarely seen ||3. Hence, 
soft quasicrystals are thought to have essentially different source of stability. 

Theoretical studies on the origin and stability of an order pattern, including periodic 
and quasiperiodic crystals, often involve minimizing a suitable free energy functional of 
the system, and comparing free energies of different candidate phases |[T4l[T^ . There¬ 
fore a systematic investigation of the stability of quasicrystals requires the availability of 
appropriate free energy functionals and accurate methods of evaluating free energies of 
quasicrystals. Several microscopic models have been studied over the years to explore the 
formation of quasicrystals arising from pair potentials with more than one microscopic 
length scales. Among them, Denton and Lowen IIT^ obtained stable colloidal quasicrys¬ 
tals within one-component system. On the other hand, phenomenological models based 
on coarse-grained free energy functionals are widely applied and particularly useful to 
the study of the stability of soft quasicrystals |[T7ll2^l^ . The earliest one can trace back 
to Alexander and McTague IITtII . They showed the possibility of stabilizing icosahedral 
quasicrystal using a Landau-type free energy functional with one order parameter. Re¬ 
lated works include Bak IITSU , Jaric Kalugin ei al 1201 and Gronlund and Mermin 1^ . 
Mermin and Troian II22II followed Alexander and McTague's theory, but introduced a sec¬ 
ond order parameter to obtain stable icosahedral quasicrystal. Swift and Hohenberg Il23ll 
added a positive-definite gradient term into the free energy functional to represent the 
effect of characteristic length scale. From the viewpoint of Fourier space, the gradient 
term is small only near the critical wave number kc — 1, thus suppressing the growth of 
any instabilities with wave numbers away from this value. Their model can be used to 
describe the supercritical instability transition from a homogeneous state to one-mode 
patterns. After that, Muller Il24ll used a set of two coupled partial differential equations; 
the pattern of a primary field is stabilized by a secondary coupled field which provides 
an effective space-dependent forcing. Two-dimensional quasicrystals, consisting of 8- 
and 12-fold orientational order have been obtained in their models. Later, Dotera 1251 ex¬ 
tended Mermin-Troian model to ABC star copolymer systems with incompressible condi¬ 
tion; several ordered two-dimensional phases were investigated, including quasicrystals 
with decagonal and dodecagonal symmetry and an Archimedean tiling pattern named 
(3.3.4.3.4). 

Beyond the study of soft quasicrystals, Lifshitz and Petrich II26II investigated qua¬ 
sicrystalline patterns arising in parametrically-excited surface waves. Although moti¬ 
vated by different physical phenomena, the free energy functional of their model is in 
a similar form with those in the preceding Landau-type models. They used only one 
order parameter describing the amplitude of the standing-wave pattern, yet an addi¬ 
tional differential term characterizing multiple-frequency forcing in the free energy func¬ 
tional. They successfully stabilized dodecagonal quasicrystalline pattern by introducing 
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multiple-frequency forcing and three-body interactions. 

Although the forementioned theoretical studies focus on different physical issues, 
they reveal the same physical nature in the formation of quasicrystals, (i) They may 
arise from the competition of multiple length scales |[T5l[T^l22ll24l - [26ll . (ii) At the same 
time the three-body interactions characterized by the cubic terms play an important role 
in stabilizing soft quasicrystals II27II . Therefore constructing an appropriate free energy 
functional to study the stability of quasicrystals should consist of more than one length 
scale and three-body interactions at least. It is also noticed that the soft quasicrystals are 
usually discovered in relatively complicated multicomponent soft-matter systems, such 
as in ABC star-shaped terpolymers IITqII , ABA'C tetrablock copolymers ||12|. Therefore a 
free energy functional with multiple coupled fields may be appropriate to study these 
complex soft-matter systems. In this work, we will make full use of those underlying 
physical nature of quasicrystals in the literature and extend the Swift-Hohenberg model 
to study the stability of soft quasicrystals. The extension is twofold: one is a positive- 
definite gradient terms with multi-mode interaction to describe competing length scales; 
the other one is that we introduce multiple order parameters which would be available 
to study the generic multicomponent systems. 

Besides a proper free energy functional, the study of thermodynamic stability requires 
accurate evaluation of free energy of various ordered phases. A number of methods have 
been developed to calculate energy of quasicrystalline pattens, including two-mode ap¬ 
proximation approach inSl26ll , using periodic structures to approximate quasicrystals 
with large unit cells Il24l - l26ll28ll29l . An alternative approach to calculate the free energy of 
quasicrystals is based on the observation that quasiperiodic lattices can be generated by 
a cut-and-project method from higher-dimensional periodic lattices ll3Ql - l32ll . It follows 
that the density and free energy of quasicrystals can be obtained using the quasiperiodic 
lattices derived from the higher-dimensional periodic structure. A approach along this 
line is the Gaussian method, in which the density profile of a quasicrystal is assumed 
to be given by a sum of Gaussian functions centered at the lattice points of a predeter¬ 
mined quasicrystalline lattice II33II . The width of the Gaussian function is treated as a 
variational parameter, which is optimized to minimize the free energy of the system. 
More recently, Jiang and Zhang 13411 proposed the projection method which embeds the 
Fourier space of quasicrystal into a higher-dimensional periodic structure, so that the 
quasiperiodic pattern can be recovered by projecting the higher-dimensional reciprocal 
lattice vectors back to the original Fourier space through a projection matrix. As a special 
case, it can also be used to investigate periodic crystals by setting the projection matrix as 
an identity matrix. From this perspective the projection method provides a unified com¬ 
putational framework for the study of periodic crystals and quasicrystals. In this work, 
we will continue to develop the projection method by optimizing the shape and size of 
higher-dimensional periodic cell which is important to explore potential metastable or¬ 
dered patterns and evaluate the free energy with higher accuracy. 

The organization of the rest of article is as follows. In Section |2l the coupled-mode 
Swift-Hohenberg model with two coupled fields is presented. We will see that no fore- 
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mentioned artificial assumption is needed in its analysis with the extension. In Section 0 
we will discuss a two-mode approximation approach and a computational framework to 
analyze the proposed model. The former mainly focus on the behavior of the model in 
the limiting regime, when strong constraints are imposed to the wave numbers; the latter 
concerns the general case. Our numerical method is based on the projection method Il34ll , 
combined with a variable cell method optimizing higher-dimensional periodic cells in 
order to deal with potentially quasicrystalline patterns and other complex phases. In 
Section IH phase diagrams for stable decagonal and dodecagonal quasicrystals in the 
coupled-mode Swift-Hohenberg model will be exhibited. We will also emphasize the 
important roles of the ratio of characteristic length scales in pattern formation. Section 0 
will be devoted to conclusion and a brief discussion. We will remark that, the numerical 
framework we develop here can well apply to a large class of more sophisticated models 
for soft quasicrystals, for example, those employing convolution-type kernels of pair¬ 
wise correlations 1^1^ . And our study of the coupled-mode Swift-Hohenberg model 
essentially promotes understanding in more general soft-matter cases. 


2 Coupled-Mode Swift-Hohenberg Model 

We consider the stability of soft quasicrystals in three-component systems, such as ABC 
star-shaped terpolymers IIIOL ABA'C tetrablock copolymers II12L in the Landau theo¬ 
retical framework. Quantities of interest are effective densities of each component, de¬ 
fined by the deviation of spatial densities of monomers from their corresponding av¬ 
erages. We denote the effective densities of components A, B, and C by and 

<I>c respectively. By assuming local incompressibility, = 0/ the free energy 

of the three-component system can be expressed in terms of two order parameters II25L 
ip—^A+^B and (p=<i>A—^B, which should appear in the Landau-type free energy 

functionals. 

Then we construct the minimal model to stabilize the periodic and quasiperiodic pat¬ 
terns. Followed by the work of Alexander and McTague lllTl . Muller ||2l|, and Dotera II25I , 
the nonlinear term in the free energy functional containing two-, three-, and four-body 
interactions is written as 

/norlmear[!/’/^] = :^ jdr + gQip^+ + t(p^+ tQ(p^+ (p^ -giXp^(p-g2^p(p^]. (2.1) 

where t, r, to, go, gi and g 2 are parameters depending on the interaction between com¬ 
ponents and the thermodynamic conditions, such as temperature. In the above free en¬ 
ergy functional, the third-order terms play an important role in stabilizing periodic and 
quasiperiodic crystals, because the cubic term is associated with triad interactions among 
wave modes. According to Dotera, when the component of C dominates the disorder- 
order transition, the term (p^ in (12.Ij) can be neglected 1125IL i.e. we set to = 0. The quartic 
terms are responsible for providing a lower bound for the free energy. Since only the 
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minimal model is constructed to study ordered patterns of interest, for simplicity, the 
quartic cross terms are omitted. 

As is mentioned in the preceding section, the formation of quasicrystals results from 
multiple characteristic length scales. In order to introduce multiple frequency forcing, 
inspired by the original work of Swift and Hohenberg 1^ , later by Muller ||24|, and 
Lifshitz and Petrich II26II . we add two positive-definite gradient terms into the free energy 
functional, acting on two order parameters respectively, 

= h{^ 

V^2J ( 2 . 2 ) 

+ J dr{rip^+gotp^ + ip^ + t(p'^ + to(p^+(p^-gitp^(p-g2q^<p^)j- 

They act as soft constraints on the magnitude of wave vectors, or the wave numbers, of 
the order parameters, with c>0 indicating their intensity. Two critical wave numbers are 
preferred, kc = l and q 1). The Fourier modes with wave numbers away from them 
will increase the free energy in the differential terms and thus are suppressed. Indeed, if 
we formally apply ParsevaTs identity, we find that 

|dr[(V2 + l),/,]'+[(V2 + ^j2)<^]'=^|dk(-|k|2 + l)2|^(k)|2+(-|k|2 + ^j2)2|^(k)|2^ 

(2.3) 

where (resp. is the Fourier transform of xj) (resp. (p) and d is the spatial dimension. 
The expression is similar in the case when xp and (p can be characterized by Fourier series 
instead of Fourier transform, simply obtained by replacing the integrations by summa¬ 
tions. As a result, in order to obtain lower energy, the Fourier wave vectors of xp (resp. (p) 
tend to be close to the circle with radius 1 (resp. q) in the Fourier space. If c ^ +oo, xp 
(resp. (^) should be strictly supported on the circle with radius 1 (resp. q). This recovers 
the forementioned assumption that the wave vectors corresponding to xp and (p should 
favor a single magnitude in much literature Il22ll25ll , and makes this model numerically 
computable. It should be noted that, there might be only one principal wave number for 
periodic structures, which implies that one of the order parameters is degenerate. For 
quasicrystals that are a class of multimode patterns and have multiple principal wave 
numbers, both of the critical wave numbers should be occupied Il2^l35l . 

3 Methods 

3.1 Two-mode approximation method in the limiting regime c — )■ +oo 

It is helpful to study the phase behavior of the coupled-mode Swift-Hohenberg model 
under the limit c — > -|-oo. As is noted above, the analysis performed for various Landau- 
type models in much literature is essentially of this type Ill4ll22ll251l26l . In the following, 
we need this type of asymptotic results to compare with those obtained under finite c's. 
Hence we discuss the two-mode approximation approach briefly. 
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For + 00 , the wave numbers of order parameter ip (resp. (p) should be strictly equal 
to 1 (resp. q). We can formally expand them as 

V’W= D *^kexp(fk-x), ^(x)= y) ^exp(/k-x). (3.1) 

|k|=l \k\=q 

(resp. K(p) is the set of spectrum points of ip (resp. (p), such that k G (resp. k G K(p) if 
and only if —k G (resp. —k G Kcp). ipy^'s (resp. ipj^'s) are complex Fourier coefficients of 
ip (resp. (p) corresponding to the wave vector k. Since xp and (p are real-valued functions, 
the coefficients should satisfy (^_k = (^k- Substituting (13.1|) into (I2.2|) , we 

immediately get the energy functional as a function of xpi^'s and 

Based on the choice of q and the rotational symmetry of interest, we will a priori select 
several candidates by determining and as in II25II . They will be specified below 
right before each numerical simulation. For simplicity, we also choose |t^k| (resp. to 
be identical for all (resp. k G K(p). Then we use the common optimization method, 

e.g., the steepest descent method, to minimize the free energy / with respect to ipi^'s and 
(^k's for each candidate. Note that t/?k's and are complex numbers, which implies 
we should not only optimize the amplitude of plane waves, but also their phase angles. 
Finally, we choose the configuration with the lowest energy as the ground state. Repeat 
the above procedure for all the parameter points of interest, we will obtain the phase 
diagram in the limit c ^ +oo. 


3.2 Direct minimization approach with the variable cell method 


In order to study the coupled-mode Swift-Hohenberg model for general c's, we propose 
the following numerical method to directly minimize the free energy functional (I2.2|) . 

We basically employ the steepest descent method in the minimization: introduce an 
auxiliary variable t and solve the following relaxation equations from some initial con¬ 
figurations ip{x,0) and (p{x,0) until the equilibrium is reached. 


^^-^^-c{V^ + lfip-{2Ttp+3goip^+^^)+2giip(p+g2(p^ 
^^-^^-c{V^+q^)^(l)-{2t(p+3to(p^+^^)+2g2lp(p+gllp^. 


These equations are solved in a pseudo-spectral approach. Take Fourier transform on 
both sides and we find 


^^-c{-\kf + lftp-{2Ttp+3goip^+4:ip^J+2gi{ip(p)''+g2{(p^J 
^ = -c{-\k\^+q^f^-{2t(p+3to(p^+4:<p^J+2g2{ip<pJ+gi{rp^J. 


(3.3) 


The differential terms can be easily evaluated by pointwise multiplication in the Fourier 
space, while for the polynomial terms, we first perform inverse Fourier transform to 
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obtain ip and (p. Then we compute them in the real space and perform Fourier transform 
to get the desired form in the Fourier space. 

We note that the pseudo-spectral method is suitable for this problem because all the 
periodic and quasiperiodic orders xp and (p will have discrete Fourier spectrum. On the 
other hand, it requires careful discretization of the Fourier space so that all the spectrum 
points can be well captured. The periodic orders are discretized in an ordinary way since 
their spectrum points locate on a periodic lattice in the Fourier space. To deal with the po¬ 
tentially quasiperiodic states, we employ the projection method ll34H . For 2 -dimensional 
quasicrystals, more than two basis vectors are used in the Fourier space to represent the 
discrete spectrum of quasicrystalline patterns. To be more precise, we need n>2 basis 
vectors in that are linearly independent in the field of rational number Q, namely, 
61 , 62 , - * • , 6 ^ G R^. Then all the discretization points in the Fourier space are encoded as 

k^'taieu aie{-N+l,---,N}. (3.4) 

i=l 

The linear independence guarantees that different (^i,- • • will give different k's. The 
choice of n and e^ closely relies on the orders of interest, especially their symmetry. In par¬ 
ticular, n = 4 for decagonal and dodecagonal quasicrystals. The basis vectors are initially 
chosen as 

6 / = ( cos sin V / = 1 ,2,3,4, (3.5) 

\ m mj 

where m = 10 (resp. m = 12 ) in the decagonal (resp. dodecagonal) case. We will discuss the 
approach of optimizing computational cell, or equivalently adjusting basis vectors later. 
Then all the computations will well resemble those in computing a periodic order using 
spectral method in n-dimensional space. Readers are referred to (Mil for more discussions 
on technical details. We only note that the computational cost is considerably large in 
such higher-dimensional computations. For example, for decagonal and dodecagonal 
quasicrystals, we use O(N^) spectrum points, where N is the number of discretization 
points on each dimension given in (|3.4|) . Hence the cost in each step of iteration in the 
steepest descent method grows as O(N^logN) since we heavily use FFT. 

It follows that Equations (|3.3|) becomes a group of ODEs for the Eourier coefficients on 
the given spectrum points. Due to numerical stiffness, semi-implicit scheme is employed 
in the iteration to avoid extremely small time step. Eor simplicity, we use fixed time step. 
Initial value of the iteration is given by approximants of potential ground states, such as 
some quasicrystalline orders of interest or some periodic orders whose Eourier spectrum 
consists of plane waves with magnitude 1 and q only. We will specify these initial values 
later right before each numerical experiment. 

As is discussed above, we choose the basis 6 / (and thus all the spectrum points of 
xp and (p) in the numerical method. However, it is often difficult to determine a priori 
where the spectrum of the ground state should locate. The choice of basis becomes an 
artificial restriction in the computations. Hence, we are going to improve the minimiza- 
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tion method with a strategy of optimizing the computational domain, which is named 
variable cell method. 

Along with (I3.3|) , we also apply the steepest descent method to optimize the basis 
vectors, i.e., 

^ = i = (3.6) 

at oei 

where A is a positive constant determined empirically to make the evolution numerically 
stable and efficient. Explicit scheme and fixed time step are employed. Since / is orig¬ 
inally not a function of the basis vectors e/s, we note that they actually come into the 
differential term of the free energy functional. To see this, we can expand ip and (p as in 
(HU). Then the free energy becomes 

k 


where c' is a positive number depending on c and the constant coming out from the 
integral, xp (resp. is a long vector whose components are all ipj^'s (resp. g(xp,^) is 
the energy contributed by the polynomial terms in (I2.2|) , which does not depend on the 
position of Fourier spectral points (and thus the basis vectors) but only their coefficients. 
One then expands k in the above formula by using (I3.4|) . It is easy to see that the free 
energy becomes a function of xp, ^ and e/s. Hence, (13.6|) is well-defined. 

There is still a special constraint in optimizing basis vectors in the quasiperiodic or¬ 
ders, when the number of basis vectors is larger than the dimension of the Fourier space: 
the basis vectors should be linearly independent in the field of rational numbers Q. Oth¬ 
erwise, the expression of the energy contributed by the polynomial terms (the second 
term in (13.7|) ) will be changed. Once the basis vectors start to evolve with (|3.6|) , it is dif¬ 
ficult to preserve or check this property. However, a weaker condition is helpful: the 
basis vectors should not be parallel with each other. This condition can be violated in 
the numerical sense. For example, when one starts minimization with a quasiperiodic 
configuration using four basis vectors under some parameters, where a two-dimensional 
periodic order is the actual ground state, the four basis vectors may tend to "mimic" the 
spectrum of the two-dimensional periodic order. Some basis vectors will thus become 
almost parallel to each other and the problem becomes ill-conditioned. To formulate this 
criterion, we let i = 1,2,...,n, and require that 


12 12 




(3.8) 


where e is a small positive number whose typical value is chosen as 0.05. In every step of 
steepest descent iteration, we will check if (13.8|) is satisfied. If not, the optimization has to 
be terminated. In computing periodic orders where the number of basis vectors is equal 
to the dimension of the Fourier space, the same criterion may also apply although it is 
not necessary. 

Based on the discussions above, we conclude the direct minimization algorithm by 
sketching the numerical recipes 
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Step 1 Given parameters c, q, i, to, t, go, gi and ^ 2 , choose the basis vectors and the initial 
configurations of xp and (p. 

Step 2 Freeze the basis vectors and minimize (I2.2|) with respect to xp and (p by (|3.3|) . The 
iteration is terminated when the right-hand side of (I3.3|) is sufficiently small. 

Step 3 Minimize (I2.2|) with respect to xp, (p and basis vectors simultaneously by (I3.3|) and 
(13.6|) . Condition (13.8|) is checked after each step of iteration. The minimization is 
terminated when the right-hand sides of (I3.3|) and (13.6|) are both sufficiently small, 
which is a good ending, or when (13.8|) is violated, which is a bad ending. In the 
latter case, a new initial configuration has to be chosen. 

Step 4 Repeat the above process with various initial configurations and find out the 
ground state with the lowest free energy. 

By repeating the above steps with different choice of parameters, we can obtain phase 
diagrams of interest. We also remark that Step 2 helps improve the efficiency of the algo¬ 
rithm, since the initial guess of basis vectors can be far away from the final configuration, 
which may lead to long iterations if we optimize order parameters and basis vectors si¬ 
multaneously from the beginning. 

In principle, the developed numerical methods can be applied to compute two- and 
three-dimensional quasiperiodic patterns. The dimension n of computational space in 
the projection method depends on the symmetries of considered quasicrystals IMlI . For 
example, for three-dimensional icosahedral quasicrystals, the projection method should 
be performed in six-dimensional space, which results in huge computational cost. Hence, 
we restrict the current study to two-dimensional quasiperiodic patterns whose point 
group symmetries can be embedded into four-dimensional periodic lattice, such as qua¬ 
sicrystalline orders with 10- or 12-fold symmetry. 


4 Numerical Results 

4.1 t — T phase diagrams in the limiting regime c — )■ +oo 

We will focus on searching for decagonal and dodecagonal quasicrystalline orders under 
the limit + 00 in this subsection. Due to the geometric features of the desired patterns, 
we will choose q = 2cos j and q — 2cos ^ for the decagonal case and dodecagonal case, 
respectively, as in 1551051 . We also let fo=0, go=0.2, gi =2.2 and g2=2.2 for the decagonal 
case, and to = 0, go — 0.8, gi = 2.2 and g 2 = 0.2 for the dodecagonal case. We leave t and t 
as free parameters forming the phase space. The above settings are the same as in I25II . 

In the decagonal case {q = 2cos{Tc/5)) and the dodecagonal case (£j = 2cos(7r/12)), we 
select K^'s and X^'s of the candidates as in FIG. [Hand FIG. Irrespectively. The black dot 
in each small figure represents the origin of the Fourier space. The blue (resp. red) dots 
surrounding the origin are end points of wave vectors in K^p (resp. Kp). These candidates 
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Figure 1: Fourier spectrum of candidates in the decagonal case: ij=2cos§, = = gi=2.2 and g2=2.2. 

The black dot in each small figure represents the origin of the Fourier space. The blue and red dots surrounding 
the origin are end points of the wave vectors in K^p and K(jj respectively. 


Figure 2: Fourier spectrum of candidates in the dodecagonal case: c] = 2cos^, tQ = 0, = Sl = 2..2 and 

g 2 = 0.2. The black, blue and red dots have the same meanings as in FIG.[TJ 

are chosen either because large numbers of triangles can be formed in the Fourier space 
using their wave vectors 1261 which can reduce the free energy value, or because their 
Fourier spectrum is similar with that of some widely observed orders in structured soft 
materials, such as lamellae phase and lamellae phase with alternating beads in block 
polymer systems 13911411 . We also remark that the first candidates in FIG. [Hand FIG. [2] 
enjoy decagonal and dodecagonal rotational symmetry respectively. 

With the two-mode approximation method discussed in Section 13.li we can obtain 
phase diagram in the t — r plane for the decagonal and dodecagonal case in the limit 
c + 00 , shown in FIG. IH Note that only one q, as specified above, is used throughout 
each individual phase diagram. The phases are named according to their real-space mor¬ 
phologies. To be more precise, we first recover the effective densities Og and Oc 
from ip and (p through 

= = (4.1) 

Then the morphology of phases are determined by plotting dominant regions of the three 
components; the dominant region of a component is defined to be the region in which its 
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effective density is the highest among the three II39II . Readers are referred to Section l42l 
for illustrations of this transformation. 



(a) Decagonal case 



(b) Dodecagonal case 


Figure 3: Phase diagrams in t — r plane under the limit c^+oo: (a) decagonal case: ^J = 2cos^, to = 0, go = 0.2, 
gl=2.2, g 2 = 2.2-, (b) dodecagonal case: ij = 2cos^, = = S’! =2.2, g 2 = 0.2. The phases are named 

after their real-space morphology. D and DD are decagonal and dodecagonal quasicrystals respectively. L+B 
is the lamellae phase with alternating beads. L 2 and /.3 are two types of lamellae phases. Hex is a hexagonal 
phase. D/s is the uniform, isotropic disordered phase. 


The phase diagram for the decagonal case (FIG. |3(a)1 ) is largely symmetric with respect 
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to the diagonal t — r. The decagonal quasicrystal, named D, is stable in a heart-shaped 
region around the critical temperatures of the two order parameters, i.e. (t,T) = (0,0). 
Some periodic orders are stable in the regions around. L+B is short for lamellae phase 
with alternating beads II39II . It is stable in four separate regions. L 3 and L 2 are two lamellae 
phases. In L3, all of the three components have their own dominant regions and they 
form lamellae phase. In L 2 , two of the components form lamellae phase, while the third 
component fully blends with them. Hex is a hexagonal phase, where one component is 
enriched in a hexagonal lattice and the other two components coexist in the remaining 
matrix. L 2 and Hex can be treat as degenerate phases, since one or more components have 
no dominant regions in the real space. They occur when t or r is sufficiently positive so 
that one order parameter has to vanish. Dzs, short for disordered phase, is a uniform and 
isotropic phase where t/? = (^ = 0 . 

In FIG. |3(b)t a small region of the dodecagonal quasicrystal (DD) is found close to 
(t,T) = (0,0). However, the region is entirely in t >0 half-plane. This implies that under the 
current parameters, the dodecagonal quasicrystalline order of (j) is completely induced by 
xp through the coupling effects. L+B, Hex and L 2 are also found in the phase diagram, with 
the same morphology as above. 

We remark that our result is qualitatively similar with Dotera's II25II in the decagonal 
case, but it is quite different in the dodecagonal case. The reasons are threefold: firstly 
we optimize all the Fourier coefficients independently rather than assuming strong ro¬ 
tational symmetry of the coefficients; secondly, we only put orders corresponding to the 
same parameter q into comparison; lastly, we do not consider more complex orders such 
as the Archimedean tiling (3.3.4.3.4). 

4.2 t — r phase diagrams by the direct minimization with the variable cell 
method 

We set c = 80 in the following simulations as an example to present the phase behavior 
under finite c's. The parameters are set as the same as before, namely, q = 2cos^, to = 0, 
go = 0 . 2 , gi = 2.2 and g 2 = 2.2 in the decagonal case, and q = 2 cos to = 0 , go = 0 . 8 , gi = 2.2 
and g 2 = 0.2 in the dodecagonal case. In order to obtain the phase diagrams, we need sev¬ 
eral initial configurations of the spectrum points. These initial configurations are selected 
based on the ground states observed in the preceding phase diagrams in the limiting case, 
i.e., FIG. |3l One reason is that we believe the phase behavior in the case of c = 80 is quali¬ 
tatively close to that in the limiting case c ^ + 00 . Also, the variable cell method can find 
the optimal configuration of spectrum points by itself, since the basis vectors are allowed 
to vary in the minimization. Hence, there is no need to choose as many initial candidates 
as in the limiting case. We represent the initial configurations again by specifying their 
Fourier spectrum, i.e., and Kcp in (13.1|) . See FIG. 0] for the decagonal case and FIG. |5] 
for the dodecagonal case. The black, blue and red dots in these figures have the same 
meanings as in FIG. [Hand FIG. [H The first candidate in FIG. |l](resp. FIG. 0 stands for 
decagonal (resp. dodecagonal) quasicrystal with perfect symmetry. We apply the projec- 
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tion method in the minimization and use four basis vectors. The other initial candidates 
are periodic orders in two dimensional space. 


Figure 4: Fourier spectrum of initial configurations in the decagonal case: ^J = 2cos^, to = 0, go = 0.2, gi=2.2 
and g 2 = 2.2. The black, blue and red dots have the same meanings as in FIG.[TJ 


Figure 5: Fourier spectrum of initial configurations in the dodecagonal case: ^J=2cos^, to = 0, = gi=2.2 

and g 2 = 0.2. The black, blue and red dots have the same meanings as in FIG.[T1 


With the numerical method proposed in Section [3^ the phase diagrams are obtained 
in both decagonal and dodecagonal case in the t — r plane with c = 80, as is shown in 
FIG.m Again, we note that only one q is used throughout each individual phase diagram 
in computing all the ordered phases. Morphology of the microphases in the both cases 
are exhibited in the FIG. [9] and FIG. [TO] respectively. As before, they are determined by 
plotting the dominant regions of components A, B and C II39II . As illustrations, we show 
in FIG. [Zl and FIG. [ 8 ] how we transform the fields of tp and (p into the effective densities 
Ob and Oc and then into the real-space morphology in the cases of decagonal and 
dodecagonal quasicrystals respectively. Readers are referred to (|4.H) and the discussions 
associated to that in Section ITT] for more details of performing this transformation. 

In the decagonal case (^J = 2 cos( 7 r/ 5 )), the new phase diagram is similar with the one 
in the limiting case shown in FIG. |3(a)| The stability region for the decagonal quasicrystal 
D slightly swells, while Hex disappears in the phase diagram. Moreover, a new transition 
phase named lamellae phase with beads at the interface, L+BI for short, is discovered 
between L+B and L 2 phases with r > 0. We also note that the L+B phase has four different 
shapes in different regions of the phase diagram, as FIGs. |9(e)l|9(h) show. 

The new phase diagram in the dodecagonal case (^J = cos( 7 r/ 12 )), i.e. FIG. | 6 (b)[ sig¬ 
nificantly differs from the one in FIG. |3(b)[ The stability region of DD shrinks while three 
new phases are observed. A new three-color lamellae phase L 3 is favored when r is suf¬ 
ficiently negative. The hexagonal phase with beads, denoted by Hex+B, appears in the 
place of Hex in FIG. |3(b)[ We note that (p = 0 in Hex when c ^ 00 , while (p^O in Hex+B for 
finite c. The morphology of Hex+B is shown in FIG. |10(c)[ The third new phase, named 
dodecagonal quasicrystalline phase with hexagonal modulation, or DD-Hmd for short. 
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(a) Decagonal case 



(b) Dodecagonal case 


Figure 6: Phase diagrams in plane for c=80: (a) decagonal case: q=2cos^, tQ=0, gQ=0.2, gi=2.2, g2=2.2] 
(b) dodecagonal case: ^j = 2cos^, to = 0, go = 0 - 8 - S’! =2.2, g 2 = 0.2. The phases are named after their real- 
space morphology. L+BI is the lamellae phase with beads at the interface, viewed as a transition phase between 
L+B and L2. Hex+B is the hexagonal phase with beads. DD-Hmd is short for dodecagonal quasicrystalline 
phase with hexagonal modulation, which can be treated as a superposition of a perfect quasiperiodic order with 
dodecagonal symmetry and a periodic order with hexagonal symmetry and also a transition phase between the 
two. Other names of the phases have the same meanings as in FIG. |3l 
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(a) D phase 


(b) Order parameter xp 


(c) Order parameter (p 
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(d) Effective density of A, (e) Effective density of B, Og 


(f) Effective density of C, 


Figure 7: Illustration of the transformation from order parameters to real-space morphology in the case of 
decagonal q uasic rystal obt ained under c = 80, q = 2cos{n/5), to = 0, = S’! =2.2, g 2 = 0.2, t = 0 and 

T = 0. FIG. |7(b)| and FIG. |7(c)| show the order parameters obtained from the dir ect m inimizati on o f the free 
energy functional (12.21) using the num erica l method pr opose d in Section [3^ FIG. r7(d)l to F IG. [7(f)] show the 
effective densities recovered through <4.1t . Then FIG. |7(a)| which is simply a copy of FIG. |9(a)|to show the 
real-space morphology of the decagonal quasicrystal, is obtained by finding out the dominant regions of the 
three components. The blue, green and red regions represent dominant regions of the component A, B and C 
respectively. 























(a) DD phase 


(b) Order parameter ip 


(c) Order parameter (p 



I 




(d) Effective density of A, (e) Effective density of B, Og 


(f) Effective density of C, 


Figure 8: Illustration of the transformation from order parameters to real-space morphology in the case of 
dodecagonal quasicrystal obtained under c = 80, q = 2cos{n/12), ^j = 2cos^, tQ = 0, = — 5’2 = 0-2, 

t=0.3 and t=-0.2. FIG. |8(b)lan d FIG. |8(c)|show the order parameters obtained from the di rect m inimiz ation 
of the free energy functional II2.2D using the numerical method pr oposed in Section [3^ FIG. |8(d)|to FI G. |8(f ) 
show the effective densities recovered through Il4.ll) . Then FIG. |8(a) which is simply a copy of FIG. 10(a 


to show the real-space morphology of the dodecagonal quasicrystal, is obtained by finding out the dominant 
regions of the three components. The blue, green and red regions represent dominant regions of the component 
A, B and C respectively. 
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(a) D phase at ^ = 0 ,t = 0 


(b) L 3 phase at t = 1 ,t = — 0.5 (c) L 2 phase at t = — 0.5,r = 1 



;»»} 




(d) L+BJ phase at t = —1.4,t = 1 (e) L+B phase at t = 0,T= —2 (f) L+B phase at t = —0.1,t = 0.7 




(g) L+B phase at t = —2,t = 0 (h) L+B phase at t = 0.7,T=—0.1 


Figure 9: All phases in the phase diagram (see FIG. 6(a)[ ) in the decagonal case {q=2cos{n/5)). The morphology 
is determined by plotting dominant regions of the three components (see FIG. [7] and FIG. [8] for illustrations). 
The blue, green and red regions represent dominant regions of the component A, B and C respectively. Note 
that L+B phases have four different real-space morphologies in different regions of the phase diagram. 
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is not a periodic order. FIG. |10(b)| shows its real-space morphology. It can be well de¬ 
scribed by a superposition of a perfect quasiperiodic order with dodecagonal symmetry 
and a periodic order with hexagonal symmetry. It arises between the stability regions 
of the perfect dodecagonal quasicrystal phase (DD) and the hexagonal phase with beads 
(Hex+B), which can be viewed as a transition phase between the two. To the best of our 
knowledge, it is the first time that such a modulated phase is reported in the numerical 
simulation. 


The phase diagram in the dodecagonal case for c = 80 clearly shows the effect of c on 
the phases behavior. In the case of finite c, Fourier wave numbers other than 1 (resp. q) 
are allowed in the spectrum of xp (resp. (p). Asa result, in Hex+B phase, p is allowed to be 
non-zero. Also, in the newly found L 3 phase at the bottom of the phase diagram FIG. | 6 (b)[ 
the magnitudes of principle wave vectors in the spectrum of xp (resp. (p) are not 1 (resp. q), 
which would be excluded in the limiting case. We also remark that if no optimization is 
performed in the basis vectors, this L 3 phase can hardly be captured. 



(c) Hex+B phase at ^ = 0.5 ,t = —0.3 (d) L+B phase at t = —0.1 ,t = 0 


Figure 10: Selected phases in the phase diagram (see FIG. | 6 (b)] ) in the dodecagonal case {q = 2cos{n/12)). 
The morphology is determined by plotting dominant regions of the three components (see FIG. [7] and FIG. [ 8 ] 
for illustrations). The blue, green and red regions represent dominant regions of the c ompo nent A, B and C 
respectively. For conciseness, the pictures for L 3 and L 2 phases are omitted. See FIG. |9(b)| and FIG. |9(c)| in 
the decagonal case for reference. 
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4.3 Phase diagram in q 

In the above discussion, we have noted that q is a system-specific parameter that de¬ 
termines the ratio of two characteristic length scales between xp and (p. It is interesting 
to study how the parameter q may influence the pattern formation. In the sequel, we fix 
c=80, to=0, go— 0.2, gi —g2—2.2, t=0.3 and r— —0.2, and apply the numerical method de¬ 
veloped in Section [3^ to find out the ground states for q G [1.00,2.10] with resolution 0.01. 
The choices of initial configurations in the computations are much similar with those in 
FIG. m and FIG.|5l The stability ranges of ground states with distinct real-space morphol¬ 
ogy are shown in FIG.fTTl Besides the patterns that have been already observed in FIG. [9] 
and FIG. [TOl the newly discovered structures are exhibited in FIG. [121 Periodic phases 
we have found before such as Hex and Hex+B, and quasicrystalline phases D and DD are 
observed once again, along with three new phases. The hexagonal phase with beads at 
the interface, Hex+BI for short, occurs as a transition between Hex phase that is stable in 
q G [1.00,1.03] and L+B at G [1.06,1.36]. Its real-space morphology is shown in FIG. |12(b)| 
Square phase with beads (SQ+B) is a periodic phase with four-fold symmetry found in 
q G [1.37,1.46]. See FIG. |12(c)| A third phase (see FIG. |12(dJ] > observed in G [1.80,1.83] is 
neither a periodic order, nor any quasicrystalline orders with rotational symmetry. We 
claim that it is a superposition of L+B phase and an incommensurate modulation, thus 
named as L+B-md. It is still not clear, however, whether such phase can be observed in 
the experiments and whether it implies another type of stable phase which we are not 
yet able to capture in the current simulations (for example, three-dimensional ordered 
phases could be favored in this range). 

L+B-md 

Hex Hex+BI 3 q+b d Hex+B / DD L 3 

I ■ H 

1 14 te ts 2 22 

q 

Figure 11: One-dimensional phase diagram as a function of q under c = 80, to = 0, go = 0.2, gi=g 2 = 2.2, t = 03 
and T=— 0.2. Phases with distinct real-space morphologies are marked with labels and different colors. The 
L+B phase is favored in all the ranges colored as black; its label is omitted for conciseness. The stability ranges 
of all the phases discovered here are also summarized Table [TJ 

To further explain the stability affected by the characteristic length scale, we define 
the critical length scale of a phase to be the length scale geometrically favored by the 
desired (i.e. standard) configuration of its wave vectors. For example, in the decago¬ 
nal and dodecagonal quasicrystalline orders, the critical length scale ^J*'s are 2cos^ and 
2cos^ respectively, see Section |4T] or Refs. (11 . Tabled] gives the critical length scale ^J*'s 
and the stability intervals of various patterns. Note that the for L+B phases is not well 
defined, since their Fourier wave vectors can be adjusted to fit a wide range of q. It is 
observed that the stability ranges of most phases in Table [T] are neighborhoods of their 
critical length scale ^J*'s. 

The ground-state energy as a function of q is plotted as the red curve in FIG. O where 
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(a) Hex phase at q = 1.00 


(b) Hex+BI phase at q = 1.05 



Figure 12: Selected phases in the phase diagram in terms of q for fixed c = 80, tQ = 0, go = 0.2, gi= g 2 = 2.2, 
t = 0.3 and T=— 0.2. The morphology is determined by plotting dominant regions of the three components 
(see FIG. [7] and FIG. [8] for illustrations). The blue, green and red regions represent dominant regions of the 
component A, B and C respectively. The stability ranges of all the phases discovered \n q ^ [1.00,2.10] are 
summarized in FIG.[TT]and Tabled 
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Table 1: Critical length scale and stability ranges in q of various phases in FIG. [TT] 


Phases 

Critical length scale 

Stability ranges in q phase diagram 

Hex 

1 

[1.00,1.03] 

Hex+BI 

N/A 

[1.04,1.05] 

SQ+B 

2cos(7r/4)«1.414 (4-fold symmetry) 

[1.37,1.46] 

D 

2 cos(7r/5) ~ 1.618 (see Ref. lH) 

[1.57,1.66] 

Hex+B 

2 cos(7r/6) ~ 1.732 (6-fold symmetry) 

[1.70,1.76] 

L+B-md 

N/A 

[1.80,1.83] 

DD 

2cos(7r/12)«1.932 (see Ref. (I) 

[1.91,1.94] 

U 

2 

[1.95,2.05] 

L+B 

N/A 

[1.06,1.36], [1.47,1.56], [1.67,1.69], 
[1.77,1.79], [1.84,1.90], [2.06,2.10] 


N/A implies that the critical length scale can not be well defined in these cases. 


the shape and size of the computational cell have been optimized by the variable cell 
method. It is clear that the energy will be significantly lowered when q is close to ^J*'s of 
the phases in Table [T] whenever is well-defined, forming "basins" on the energy curve. 
As a comparison, we also plot the energy obtained without optimizing computational 
domain as the blue curve in the same figure. We can see that the red curve is always below 
the blue one, and it makes considerable improvement in evaluating energy in the basins 
compared to the blue curve. This implies that the variable cell method can better capture 
phase behavior of the model than a simple minimization. The reason is that, when q is 
slightly inconsistent with the of a particular phase, the energy value of the phase can 
be lowered by slightly adjusting the basis vectors, especially their lengths. When q is far 
away from however, chances become small such that the phase corresponding to 
can be distorted to fit q. 

5 Conclusion and Discussion 

In this article, we develop a coupled-mode Swift-Hohenberg model with two order pa¬ 
rameters to study the stability of soft quasicrystals in three-component systems. Two 
characteristic length-scales are added into the original Swift-Hohenberg free energy func¬ 
tional through the positive-definite gradient terms acting as penalizing factors on Fourier 
wave numbers. In order to investigate ordered structures, a computational framework is 
also proposed combining the projection method and the variable cell approach. It enables 
accurate calculations of the free energy of ordered patterns, leading to phase diagrams 
with more accurate phase boundaries. Moreover, it shows the potential to discover more 
complex phases, such as incommensurate modulation structures. The generalized model 
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Figure 13: Curves of the ground-state energy as functions of q for c=80, tQ=0, gQ=0.2, gi=g2=2.2, t=0.3 and 
T=— 0.2. The blue curve represents the energy obtained without optimizing computational cell, while the red 
one is the improved energy curve after optimizing computational cell. One can see the strategy of optimizing 
computational cell makes improvement in energy in the basins. 


as well as the numerical framework provide a powerful tool for investigating the phase 
behavior of ordered patterns on the mean-field level. 

Employing this tool, we study the phase behavior of the coupled-mode model under 
the soft constraints on wave numbers (i.e. c = 80), and compare it with that in the limiting 
case c +00 (i.e. infinitely stiff constraints), which is obtained by the two-mode approx¬ 
imation method. The value of q's are selected to form the decagonal (2cos(7r/5)) and 
dodecagonal (2cos(7r/12)) symmetries, respectively. Due to high computational cost, 
we only considered the two-dimensional periodic structures and quasicrystals whose 
Fourier spectrum can be embedded into a periodic lattice in the four-dimensional space. 
A number of two-dimensional ordered structures emerge in our calculations, including 
periodic crystals, 10-, 12-fold symmetric quasicrystals and modulated structures. By com¬ 
paring the free energies of these candidates, phase diagrams are obtained in the t — r 
plane in both cases of c = 80 and c ^ +oo. It is predicted that decagonal quasicrystal 
(^J=2cos(7r/5)) and dodecagonal quasicrystal {q=2cos{n/12)) can become ground states 
for both c's. The differences of the phase diagrams in these two cases are also emphasized 
to show the important role of c and the optimization of the computational domain. In 
particular, we remark that the classic two-mode approximation analysis a priori assumes 
the order parameters favor specific length scales k=l and k = q, and is not sufficient to 
capture the true phase behavior of models of this type. 

Furthermore, the role of the ratio of two characteristic length-scales, q, on the pattern 
formation is also studied. A one-dimensional phase diagram as a function of q with other 
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system parameters fixed is obtained, as is shown in FIG. [TTl From the phase diagram, 
we can see that the pattern formation depends on the critical length scale which is 
determined by the standard configurations of wave vectors of the corresponding ordered 
patterns. For most ordered patterns, their stable regions are simply neighborhoods of 
their corresponding critical length scales. A modulated structure is discovered for q G 
[1.80,1.83]. It could be an indication of new phases in experiments or shows the limitation 
of the current study; for example, a three-dimensional structure might be favored over 
all the two-dimensional structures in this range, which we have not studied. It should be 
emphasized that when varying the value of q, the variable cell method plays an important 
role of accurately calculate the free energy. 

These results provide a good understanding of the rich phase behavior on the mean- 
field level in the coupled-mode Swift-Hohenberg model for three-component systems. 
It confirms that the existence of two characteristic length scales and sufficiently strong 
three-body interactions can account for the stability of soft quasicrystals II27II . In the cur¬ 
rent work, the 10-, 12-fold symmetric quasicrystals, and some modulated structures can 
become stable phases. To this date, in the most experiments of soft matters, quasicrystals 
have been observed in dodecagonal symmetry, whereas decagonal quasicrystal have not 
been reported. However, it is still worth expecting the decagonal symmetric quasicrys¬ 
tals, and even modulated structures in the multi-component soft-matter systems. 

The numerical framework we developed, combining the projection method with the 
variable cell method, is suitable for studying periodic and quasiperiodic orders in the 
current model. It also works for a larger class of mean-field energy functionals. For 
example, if the differential terms in (I2.2|) are replaced by convolution-type interaction 
kernels Il36ll37ll 



(5.1) 



our method works as well. Indeed, by ParsevaTs identity, the convolution can be rewrit¬ 
ten as 



where C > 0 is a constant coming out from the Fourier transform. It is essentially in the 
same form as (I2.3|) . which can be dealt with in exactly the same way. 

A potential drawback of our numerical method is that, although it allows to perform 
minimization with larger degrees of freedom, it is not of high accuracy, especially when 
c and N are small. It is well-known that the ordinary spectral method enjoys spectral 
accuracy in dealing with periodic problems in general, but it might be not the case in 
our quasiperiodic settings. The part of energy, for which the wave vectors in (I3.4|) do 
not account, decreases slowly as N grows. In other words, there is still a considerable 
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amount of energy contributed by the wave vectors with high indices, which are not nec¬ 
essarily of high frequencies. In the current work, we set c = 80 and use N=16 bases in the 
direction of each basis vector. It is sufficient to obtain energy of the ground states with 
four- to five-digit precision and thus can be used in constructing the phase diagrams. 
However, when c is smaller, we need larger N, which is not always affordable for four- 
or higher-dimensional computations in the projection method. Recall that the computa¬ 
tional cost in one step of steepest descent iteration grows as O(N^logN) in computing 
decagonal and dodecagonal quasicrystals. The problem of low accuracy can be model- 
specific, since it does not occur in our earlier work on the Lifshitz-Petrich model Il35]l . In 
fact, the penalizing effect of the differential terms in the Lifshitz-Petrich model is much 
stronger than that in (I2.2|) . The former one is proportional to |k|^ as |k| ^ +oo, while the 
latter one grows as |k|^. This implies that the Fourier wave numbers far away from the 
prescribed length scales (i.e., 1 and q) are suppressed more strongly in the Lifshitz-Petrich 
model than in the current model. As a result, their contribution to the energy becomes 
less significant in the Lifshitz-Petrich model. Nevertheless, it is still unknown how to 
achieve higher accuracy in the quasiperiodic setting in the current setting. 
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